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Abstract 

A general set of methods is presented for calculating chemical potentials in solid and 
liquid mixtures using ab initio techniques based on density functional theory (DFT). The 
methods are designed to give an ab initio approach to treating chemical equilibrium between 
coexisting solid and liquid solutions, and particularly the partitioning ratio of solutes be- 
tween such solutions. For the liquid phase, the methods are based on the general technique 
of thermodynamic integration, applied to calculate the change of free energy associated 
with the continuous interconversion of solvent and solute atoms, the required thermal av- 
erages being computed by DFT molecular dynamics simulation. For the solid phase, free 
energies and hence chemical potentials are obtained using DFT calculation of vibrational 
frequencies of systems containing substitutional solute atoms, with anharmonic contribu- 
tions calculated, where needed, by thermodynamic integration. The practical use of the 
methods is illustrated by applying them to study chemical equilibrium between the outer 
liquid and inner solid parts of the Earth's core, modelled as solutions of S, Si and O in 
Fe. The calculations place strong constraints on the chemical composition of the core, and 
allow an estimate of the temperature at the inncr-core/outer-core boundary. 

1 Introduction 

We present here a set of techniques that allow the ab initio calculation of chemical potentials 
in solid and liquid solutions, and hence the ab initio treatment of chemical equilibrium between 
solid and liquid phases. There are many areas of chemical physics where such techniques might 
be important, but we believe they have a particular role to play in studying the partitioning of 
impurities between different phases under extreme conditions, where experiments are difficult 
or impossible. As an illustration of the power of the techniques, we will describe how we have 
applied them to study chemical equilibrium between the solid and liquid parts of the Earth's 
core. 

The techniques to be presented form a natural sequel to recent developments in the ab initio 
thermodynamics of condensed matter based on electronic density-functional theory (DFT) Q. 
For many years, DFT has been used to calculate the phonon spectra of perfect crystals Q, 
and it is only a short step from that to the calculation of free energies and other thermody- 
namic quantities in the harmonic approximation. There have already been several reports of 
DFT calculations of high-temperature crystal thermodynamics |5|, ^, 0, H], including 



solid-solid phase equilibria ||T^, 11, 12, ^] using this approach. For liquids, ab initio ther- 



modynamics first became possible with the Car-Parrinello technique [14| of DFT molecular 
dynamics simulation, which immediately gave a way to calculate such quantities as pressure, 
internal energy and temperature of a liquid in thermal equlibrium. The first DFT treatment 



of solid-liquid equilibrium was achieved by Sugino and Car [^], who used thermodynamic in- 
tegration to compute the free energies of solid and liquid Si, and hence the melting properties 
of the material. Closely related is the work of de Wijs et al. |16] on the melting of Al. We have 
recently reported DFT calculations of the free energies and melting curves of Fe @, mi and 
Al jl^ over a wide range of pressures; Jesson and Madden have recently presented ab initio 
calculations of the zero-pressure melting properties of Al using their 'orbital free' approach. 
The work of Smargiassi and Car |21] and Smargiassi and Madden [ p^ on the free energy of 
formation of defects in crystals is also relevant to the ideas to be presented here. 

Thermodynamic integration has been the key to calculating the ah initio free energy of 
liquids and anharmonic solids, and hence to the treatment of solid-liquid equilibrium. It 
provides a means of computing the difference of free energy between the ah initio system 
and a reference model system whose free energy is known. We will show that it is also the 
key to calculating ah initio chemical potentials of liquids and anharmonic solids, but here 
it is used in a rather different way. The chemical potential of a species is the free energy 
change when an atom of that species is added to the system. The difference of chemical 
potentials of two species is therefore the free energy change when an atom of one species 
is replaced by an atom of the other, or equivalently when one atom is transmuted into the 
other. The role of thermodynamic integration here is to provide a way of calculating the 
free energy change associated with such transmutations, and we shall show how this can be 
accomplished in practical DFT simulations. This general approach is closely related to ideas 
that have been used for a long time in classical simulation. A recent example of classical 
thermodynamic integration with molecular transmutation to calculate solvation free energies 
in aqueous solution can be found in Ref. which gives references to earlier literature. 

Although the techniques we shall present are fairly general, we do impose two restrictions 
at present: First, the theoretical framework is developed for the case of a two-component 
mixture; second, one of the components is present at low, but not very low, concentration, in 
a sense to be clarified in Sec. The situation envisaged therefore consists of fairly dilute solid 
and liquid solutions in coexistence. 

There are vast numbers of problems both in the chemical industry and in the natural 
world that depend on the partitioning of chemical components between coexisting phases, 
and the ability to calculate chemical potentials ah initio should make it possible to address 
some of these problems in a new way. Our original incentive for developing these techniques 
was the desire to understand better the chemistry of the Earth's core, and this is a good 
example of a problem where ah initio calculations can supply information that is difficult to 
obtain experimentally because of the extreme conditions of temperature (T ~ 6000 K) and 
pressure {p ~ 330 GPa). The core is composed mainly of Fe, and comprises an outer liquid 
part and an inner solid part [24|. The density of the outer core is ~ 6 % too low to be pure 
Fe 1 24, 25, 26, and cosmochemical and geochemical arguments show that the main light 

impurities are probably S, O and Si [^]. The inner core has grown over geological time by 
crystallisation from the outer core, and the partitioning of impurities between liquid and solid 
is crucial for understanding the evolution and contemporary dynamics of the core. The size of 
the density discontinuity (ca. 4.5 %) [29, 30| at the inner-core/outer-core boundary can only 
be interpreted if one understands this partitioning, and also provides a constraint on possible 
chemical compositions. We shall show how our ah initio techniques for calculating chemical 
potentials shed completely new light on this problem. Brief reports of these calculations have 
already appeared |27, 31, 32 1. 

In developing the theoretical basis of our techniques, we define our technical aims in Sec. |2| 
by summarising the standard thermodynamic relations describing phase equilibrium. The dif- 
ference of chemical potentials of solute and solvent atoms, and the free energy change associated 
with the transmutation of solvent into solute are discussed in Sec. ^. In Sec. ^, we then develop 
the ah initio techniques themselves. We shall explain how thermodynamic integration can be 
used to perform the solvent-solute transmutation so as to obtain the difference of chemical 



potentials in the liquid; we also describe the techniques for calculating chemical potentials in 
solid solutions, both in the harmonic approximation and using thermodynamic integration to 
handle anharmonicity. Sec. ^presents our results for the case of S, O and Si dissolved in solid 
and liquid Fe under Earth's core conditions, and summarizes the implications of the results 
for the partitioning of these impurities between the inner and outer core and the chemical 
composition of the core. Discussion and conclusions are given in Sec. |6|. 



2 Chemical equilibrium: thermodynamics 

Our task in this Section is to identify the thermodynamic quantities that will need to be 
calculated ab initio. Chemical equilibrium between two multicomponent phases is characterized 
by equality of the chemical potentials of each component in the two phases. For a two- 
component solution consisting of solute X dissolved in solvent A, equilibrium between solid 
and liquid phases requires that: 

f^xiP, Tm, Cx) = Mx(P> Tra, C^) , (1) 

T^, 4) = fi%{p, T^, 4) , (2) 

where and /^a are the chemical potentials of solute and solvent, p is the pressure and cx 
is the mole fraction of solute, with superscripts / and s for liquid and solid respectively; 
is the melting temperature, i.e. the temperature at which the liquid and solid solutions are 
in equilibrium, which depends on the impurity mole fractions. The two equations impose two 
relations between c^, and at the given p. In the low-concentration limit cx 0, /ix 
diverges logarithmically, and it is useful to write: 

^x(p, T, cx) = k-BT In cx + /ix(p, T, cx) , (3) 

where ilx{p, T, cx) is well behaved for all cx- In an ideal solution, /2x is independent of cx, but 
in reality the interaction between solute atoms causes it to vary with cx- Combining Eqs. (||) 
(^), we obtain: 

4/4 = exp[(/i^ - fl'x)/kBTrr,] , (4) 

so that the ratio of the mole fractions and 4 solid and liquid is determined by the 
thermodynamic quantities and fl^. The melting temperature entering Eq. (^) differs 
from the melting temperature of the pure solvent, and may be regarded as determined by 
Eq. (D. 

We now develop a practical way of solving Eqs. (|^) and (Q). We are interested in the case 
of moderately low cx, but we wish to take account of the variation of jlx with cx to lowest 
order. We therefore expand fix as: 

flx{p, T, cx) = /ix(p, T) + Xx{p, T)cx + 0{cx) , (5) 

and we shall systematically neglect the terms O(cj^). Since it will be important later, we stress 
here that this represents the concentration dependence of fix at constant pressure. Eq. (Q) 
then becomes: 

c|,/c^ = exp - fi^^ + A^^c^ - A5,cf,)/A:Br„] . (6) 

To obtain an equation for Tm, we need the corresponding expansion for /xa. We use the 
Gibbs-Duhem equation: 

CAdfiA + cxdfj.x = , (7) 

which gives: 

fiAip,T,cx) = ^^l{p,T) + {kBT + Xx{p,T)) ln{l - cx) + \xip,T)cx + 0{ci) , (8) 



where /i^ is the chemical potential of the pure solvent, and we have used the fact that ca = 
1 — cx- To linear order in cx, this gives: 



^a(p, T, cx) = fJ-AiP, T) - kBTcx + 0(cx) • (9) 

We apply this in Eq. (|2|) by expanding ^\{p,T) to linear order in the difference — 
between Tm and the melting temperature of pure solvent. This yields: 
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-A:BT„4 + /.?f(p,rO) + (r^-r°)(^) . (lo) 

V W T=TO 

Since /u5f(p, T^) = fi^{p,T^), we can rewrite this equation as: 

ksTmA + {Tm - T^)sf = k^TmC^ + {Tm - T^)sf, (11) 

where s% = - {dn\/dT)j,^ _rpQ is the entropy per atom of pure solvent at the melting temper- 
ature. The shift of melting temperature due to the presence of the solute is then: 



(^^--?^-^) = ^(4-ck), (12) 



kB^n 

„0 

where As^ = — s^f is the entropy of fusion of pure solvent. Eqs. (|6|) and (|l2|) must be 
solved self-consistently. 

We see from Eqs. (^) and (|l^) that the main thermodynamic quantities to be calculated 
ab initio are fJ^ and Ax for the solid and liquid solution. We also require ab initio values for 
the melting temperature and entropy of fusion of pure solvent. In addition, we shall find it 
necessary to obtain ab initio values of the partial molar volumes of the solute, which will be 
discussed later. We next turn to the statistical-mechanical considerations needed to develop a 
strategy for calculating ^ujr and Ax. 



3 Interconversion of solvent and solute: statistical mechanics 

The chemical potential of a solute X in solid or liquid solvent A is the change of Gibbs free 
energy when one atom of X is added to the system at constant pressure and temperature. In our 
practical ab initio calculations, we work at constant volume rather than constant pressure, so 
we prefer the equivalent statement that the solute chemical potential is the change of Helmholtz 
free energy when the solute atom is added at constant volume and temperature. However, it is 
impractical to add solute atoms to a system in an ab-initio simulation. It is more convenient to 
convert solvent into solute, which means working with the difference /Ux — /"a of the chemical 
potentials. The chemical potential /ix and hence jlx can then be obtained from /ix — fJ-A by 
making use of the Gibbs free energies of pure solid and liquid A, calculated separately. We 
have sketched this technique briefly in our previous papers p^, P, H, and we now describe 
it in more detail. 

The Helmholtz free energy of a system containing A^a solvent atoms and Ax solute atoms 

is: 

F{Na,Nx) = -^BTln| ^^^^^^^\^^^^^^ |^dRexp[-/3C/(AA, Ax;R)]| , (13) 

where /3 = l/k^T, and Ax and Ax are the thermal wavelengths of A and X, given by 
Aa = /i/(27rMA/cBr)^/2, with Ma the atomic mass of A, and similarly for Ax. The quan- 
tity U{Na, Nx'j'R-) is the total energy function of the system of Aa solvent and Ax solute 



atoms, which depends on the positions of all the atoms, indicated by R, and Jy dR indicates 
integration over the whole configuration space of the system contained in volume V. 

The difference of chemical potentials hxa = fJ'X — l^-A is equal to the change of F when a 
single atom of A is converted into X, and is given by: 

PLXA = F{Na-1,Nx + 1)-F{Na,Nx) 

= -kBTln{NA/Nx) - kBTln (Ai/A| 



We express this as: 



, ^, / Jy dRexp [-PU{Na -l,Nx + 1; R)] ] 
^ dRexp [-/3t/(iVA,iVx;R)] j ' ^ ' 



fJ-XA = kBTln— ^ h 3A;bT In ^ + m(cx) , (15) 

1 - cx Aa 



where we define: 



' \ /v dRexp [-;3(7(Wa,A'x;R)1 J 

The intensive quantity m{cx) depends only on pressure, temperature and concentration (we 
write it as m{p,T,cx)), or alternatively on volume, temperature and concentration (we then 
write it as 'm{v, T, cx), where v is the mean atomic volume V/{Na + A'x))- Expanding Eq. (|l^) 
to linear order in cx, we have: 

Ax 

fJ-XA = fceTlncx + kBTcx + S/cbTIu h m(cx). (17) 

Aa 

Compare now with Eqs. (^), (^) and (|9|): 

tJ'XA = ^bT In CX + ^Jix + Xxcx - /^a + ^b^cx , (18) 

and we have: 

^'(cx) +3A;Brin— = /u]^ + Axcx- (19) 
Aa 

If values are available for m{p,T,cx) at different values of cx for a given pressure p, we can 
obtain the quantities /Xjr^ = //^ — /^a — 3/!;BTln(Ax/AA) and Ax for that pressure. We then 
need the pure-solvent chemical potentials /.i^ for liquid and solid, whose ab initio calculation 
has been described in detail elsewhere ||9|, |l^, p^ . 

We conclude this Section by rewriting c|r/cx from Eq. (^) in terms of /^xa' and fi^ for 
liquid and solid: 

_ exp [[n^xA - MxA + ^xcx - ^x4) /^bT^] 
'^^'^ " exp [(/.Of - /.oj) /kBT^] • 

If the concentrations are small enough for the difference between and to be negligible, 
then = ij,^, and the denominator is unity; but in general deviations of the denominator 
from unity should be included. 



4 Ab initio chemical potentials 
4.1 The liquid solution 

For the liquid, we calculate the quantity m{cx) of Eq. ( |l6|) by a form of 'thermodynamic 
integration'. We first outline a simple way of doing this that is correct in principle, but suffers 



from practical problems; we then show how the method can be modified to give a practical 
procedure. 

Thermodynamic integration [^] is a general technique for computing the Helmholtz free 
energy difference Fi — Fq of two systems containing the same number N of atoms, but having 
different total energy functions Ui(R) and C/o(I^)- The difference Fi — Fq is the reversible work 
done on continuously switching the total energy function from Uq to Ui at constant volume, 
which is given by 

Fi-Fo= f dX {Ui - Uo)x , (21) 







where the average ( • ) is calculated in thermal equilibrium for the system governed by the 
'hybrid' energy function U\ = {1 — A)C/o + AC/i. This is a well established technique for the ab 
initio calculation of liquid-state free energies which was used in our recent ab initio 

investigation |l^ of the high pressure melting curve of Fe. 

In order to compute rn-(cx), we could in principle choose Uq to be the total energy for the 
system of A^a atoms of solvent and A'x of solute, and Ui to be the same for A^a — 1 atoms of 
A and Nx + 1 of X. We evaluate {Ui — Uo)x by performing ab initio molecular dynamics with 
time evolution generated by Ux, and taking the time average of Ui — Uq. This is repeated for 
several values of A, and the integration over A is done numerically. This type of 'alchemical' 
transmutation of A into X obviously does not correspond to a real- world process, but in terms 
of ab initio statistical mechanics is a perfectly rigorous way of obtaining the quantity m(cx). 
It demands an unusual kind of simulation: For the atom positions ri, . . .ttv at each instant 
of time, we have to perform two independent ah initio calculations, one for each chemical 
composition. As well as Uq and Ui for the given positions, we calculate two sets of ah initio 
forces Foi = — VjC/o and Fij = — VjC/i, and the linear combinations Yxi = (1 — A)Foi -|- AFij 
are used to generate the time evolution. 

The major problem with this scheme is one of statistics. The thermal average {■)x is 
evaluated as a time average, but since only a single atom is transmuted the scheme does 
not benefit from averaging over atoms. The efficiency of the averaging can be considerably 
improved if one is prepared to transmute several atoms simultaneously. If we do this, then 
instead of obtaining m(cx) at a given mole fraction cx, we obtain an integral of m(cx) over a 
range of cx values. The information we need can still be extracted, as we now describe. 

Consider the change of Helmholtz free energy when we start from A^ atoms of pure solvent 
and transmute A'x of them into solute atoms at constant volume and temperature. This 
can clearly be calculated by thermodynamic integration using the procedure outlined above. 
Denoting this change of free energy by W{N,Nx)., we can express it as: 

{ Jv exp [-(3U{N, 0; R)] J 

We then have: ^ 

W{N,Nx)= f dX{Ui-Uo)x, (23) 
Jo 

with Ui (R) = U{N - Nx, Nx; R) and Uo{R) = U{N, 0; R). Our procedure will be to calculate 
T^(A^, A'x) at several values of Nx/N = cx at a chosen volume, and then fit the results in the 
following way: 

W{N, Nx)/Nx = a + bcx. (24) 

The information needed can now be extracted by noting that for a given mean atomic volume 
V, the quantity m{v,T,cx) is: 

m{v, T, cx) = {dW/dNx)v,T = a + 2bcx ■ (25) 

It follows immediately that: 

/it A = lim m{v,T,cx) = a . (26) 

cx-^O 



To obtain Ax from the coefHcient b, we note that Ax = liinc^^o{dm{p,T, cx) / dcx)p and 2b = 
limc^^Q{dm{v,T,cx/dcx_)v- The fact that one derivative is isobaric and the other isochoric is 
significant. The quantity Ax that we seek describes the isobaric concentration dependence of 
solute chemical potential. But since our ab initio calculations are done at fixed volume, the 
immediately available quantity b is an isochoric derivative. 

The relation between the constant-pressure and constant-volume derivatives of m is exam- 
ined in the Appendix, where we show that: 

{dm/dcx)p = {dm/dcx)v - nBrivx - v^f , (27) 

where Bt is the isothermal bulk modulus, and -ux and are the partial atomic volumes of 
solute and solvent. We conclude that: 

Ax = 26 - nBrivx - va? ■ (28) 

Here, the quantities Bt^ vx and f a can be evaluated at infinite dilution. The calculation of Bt 
and v\ involves only ab initio m.d. simulations on the pure solvent, and presents no problem. 
We return below (Sec. ^ ) to the ab initio calculation of vx- With this, we have a complete 
procedure for determining and Ax- 



4.2 The solid solution 

If anharmonic effects are negligible, then the free energy of the solid can be obtained from ab 
initio phonon frequencies, so that thermodynamic integration is not needed, and no statistical 
averaging is involved in the ab initio calculations. There is then nothing to prevent us from 
calculating m(cx) directly from the free energy change when solvent atoms are replaced by 
solute atoms. We assume for the moment that this is adequate, and return later to the question 
of anharmonic effects. 

We start by considering the zero-concentration limit of m(cx), namely /^xA' which is the 
non-configurational free energy change when an atom in the perfect A crystal is replaced by 
an X atom. This can be written as: 

t tperf , tharm /r,n\ 

/^XA = /^XA +/^XA ' (29) 

where fJ-x'A^ is the free energy change for the perfect non-vibrating crystal, and /ixA^^™ 
harmonic vibrational part - we refer to 'free energy' to allow for the possibility of 

thermal electronic excitations, which are important in high-temperature Fe The calculation 
of ^xa'^^ is straightforward, and involves only the difference of ab initio (free) energies of the 
static fully relaxed crystal containing a single substitutional X atom and the static perfect 
crystal, the two systems having the same volume. 

In the high-temperature limit, where T is well above the Debye temperature, /^x'a™ can 
be written as: 

/^XA™ = ^bT^ In (ujjujn) , (30) 

n 

where oj'^ and cOn are the harmonic frequencies of the normal modes of the impure and pure 
crystals, and the sum goes over all modes. The frequencies are calculated ab initio, and we use 



the 'small-displacement' method described in detail elsewhere [^, 34, 35|. This involves DFT 
calculations of the force on every atom in the system induced by displacement of a single atom, 
and this has to be done for all symmetry inequivalent atoms and displacements. To obtain 
Ax in the harmonic approximation, we must include the effect of interactions between solute 
atoms. The key to this is to note that the calculation of the partition function, i.e. the integral 



over configuration space of Eq. (12), can be broken into (a) a sum of distinct configurations, i.e. 



assignments of solute atoms to lattice sites, and (b) an integral over vibrational displacements 



of the atoms away from their relaxed equilibrium positions for each such configuration. This 
means that the statistical mechanics of the solid solution maps exactly onto that of a lattice 
gas, and the free energy of the solid solution is: 

FiNA,Nx) = -kBT In (^e-^''-'^ , (31) 

where the sum goes over all distinct configurations 7, and is the non-configurational free 
energy of the system for each 7. 

It is convenient to relate F{Nj^, Nx) to the free energy i<A of the pure A crystal having the 
same number of lattice sites. The difference AF{N,Nx) = F{N\, Nx) — F\ is the change of 
free energy when A'x atoms of A in the pure crystal are transmuted into X atoms. We have: 



AF(iVA,A^x) = -fcBTln(^e- 

V 7 



(32) 



Now in the limit cx 0, we can neglect the interactions between X atoms, and we get: 

- Fa ^ Nxfilj, . (33) 

At higher concentrations, we need to include the free energy of interaction between pairs of X 
atoms, and we write: 

- Fa ~ Nxfi^A + lY.'t>^n, (34) 

where (j)mn is the non-configurational free energy change when a pair of X atoms are brought 
from widely separated sites in the otherwise perfect crystal to the sites m and n. We then 
have: 



AF{N, Nx) = iVx/iL - ^BFln ^ exp 



m^n j 



(35) 



In the later practical calculations, we approximate by setting (^mn equal to zero except when 
m and n are nearest-neighbor lattice sites, the interaction free energy being then called simply 

It is now an exercise in the statistical mechanics of lattice gases to show that the leading 
order in cx'- 

AF(N,Nx)=NxP^XK + iVA:Br[cxlncx + (l-cx)lncx] + 



+ -NkBTci^z[\-e-^0) , (36) 

where z is the coordination number of the lattice. The derivative dAF{N, Nx)/dNx gives us 
f^x — fJ-A, from which we straightforwardly extract Ax, which is given by: 

Ax = A^bTz (1 - e-^<^) . (37) 

As in the case of the liquid, this formula should be corrected from constant volume to constant 
pressure, so that the correct formula is: 

Ax = kBTz (1 - e-"*) - nBrivx - va? , (38) 

with Bx, vx and va the isothermal bulk modulus and partial molar volumes in the dilute 
solid solution. The calculation of Bt and va presents no problems; we return to the ab initio 
calculation of vx in Sec. [4.3| . 

In addition to using this analytic derivation to obtain Ax, we have also performed Monte 
Carlo calculations on the lattice gas to obtain numerical values of AF{N, Nx)- These serve 



both to confirm the correctness of the analytic result in the region of low cx and also to assess 
deviations from the linear dependence of jlx on cx as cx increases. 

The remaining task is to calculate the nearest-neighbor interaction free energy (p. This 
follows exactly the scheme for calculating /^x^^, where now (p is the non-configurational free- 
energy change when two neighboring atoms in the perfect crystal A are replaced by X atoms, 
minus twice /Uj^^- This can be written as: 

where (j)^'^^^ is the (free) energy change for the perfect non-vibrating crystal, and 0^^''™ is the 
harmonic vibrational part. The static part (j)^^^^ is obtained from the difference of ab initio free 
energies of the relaxed equilibrium system containing neighboring X atoms and the perfect pure 
A lattice. We obtain t/)'^'^™ fr om the harmonic vibrational frequencies of the system containing 
the neighboring X atoms by a formula exactly analogous to Eq. (pO|). 

We now return very briefly to the question of anharmonicity. In many cases, very high 
precision may not be needed for the chemical potentials, so that anharmonic corrections to 
/ix can be neglected. But in one case that will be important later, that of substitutional O 
in Fe, we know that anharmonic effects are large. The techniques we have used to treat them 
are described in detail elsewhere [31|. The strategy is based on thermodynamic integration 



between reference models representing both the pure Fe and the impure Fe/X systems, followed 
by further thermodynamic integrations between the ab initio and reference systems. 



4.3 Partial molar volumes in the liquid and solid solutions 

The partial molar volume vx of solute or va of solvent is the change of volume of the system 
when one atom of X or A is added at constant pressure and temperature. The volumes are 
related to the chemical potentials by: 

vx = (9/Ux/9p)t,cx > ^^A = ((?/UA/9p)r,cx • (40) 

We note that the total volume of the system is given hy V = Nxvx+Na.va- As for the chemical 
potentials, we find it easier to consider the interconversion of solvent and solute, and to work 
with the difference vxA = vx — va- The liquid is treated by ab initio m.d., in which the pressure 
for a given volume is calculated during the simulation. (In our practical calculations, we work 
at constant V.) The straightforward way of obtaining the dilute limit of vx is therefore to 
calculate the change of pressure Ap resulting from the replacement of a chosen number A'x of 
atoms in the pure solvent by X. The pressure change per atom dp = Ap/Nx then gives us vxA 
by the relation vxA = V5p/ Ex- 
it is clearly possible to follow the same route for the solid. However, if the solid is treated 
by harmonic frequency calculations with or without thermodynamic integration for the an- 
harmonic contribution, then the partial molar volumes must be obtained from the chemical 
potentials via Eq. (|40|). This requires calculation of /ix at different volumes followed by nu- 
merical differentiation. 



5 Illustration: chemical equilibrium in the Earth's core 

In applying the techniques to study chemical phase equilibrium between the Earth's inner 
and outer core, our aim is to show how they can yield important new information about 
the chemical composition and temperature of the core, both of which are controversial. Our 
strategy exploits the fact that that the density as a function of depth in the core is accurately 
known from seismic measurements |^^; in particular, it is quite well established that there is 
a density discontinuity of 4.5 ± 0.5 % across the inner-core/outer-core boundary (ICB) [ p9[ . 
Recent ab initio studies of the melting properties of pure Fe concur in giving a volume of 



fusion of ~ 1.8 % flSl , p^ , which is clearly much smaller. This means that there must be a 
substantial partitioning of light solute elements from solid to liquid to account for the large 
observed discontinuity. 

This discontinuity can be studied with our methods. If we suppose initially that the core 
is a binary solution of Fe with one of the leading impurity condidates S, Si or O p^ , then the 
solute mole fraction in the liquid core can be fixed by requiring that the density reproduce the 
seismically observed density. Calculation of the chemical potentials in the liquid and solid 
then gives us the mole fraction in the solid, from which we can deduce the solid density, and 
hence the density discontinuity. Agreement or disagreement with the known discontinuity puts 
a constraint on the composition. At the same time, the shift of melting temperature given by 
Eq. ([T^ ) gives us information about the temperature at the ICB. 

In the following, we first summarize the general techniques used in all the calculations 
(Sec. 5^). We then describe separately the calculations on the liquid and solid alloys (Sec. |5.2| 



and^^ respectively), presenting results for the chemical potentials and partial molar volumes. 



5.4, we then combine the results with seismic data to obtain constraints on the chemical 



In Sec. 

composition and temperature of the Earth's core. 



5.1 General techniques 

Our ab initio calculations are based on the well established DFT methods used in virtually 
all ab initio investigations of solid and liquid Fe 17, 18, ^ 



39], including our own 



previous work on pure Fe and its solid and liquid alloys with S and O |27, 31, 40, 41]. We 
employ the generalized gradient approximation for exchange-correlation energy, as formulated 
by Perdew et al. ]42], which is known to give very accurate results for the low-pressure elastic, 
vibrational and magnetic properties of body-centred cubic (b.c.c.) Fe, the b.c.c. — > h.c.p. 
transition pressure, and the pressure-volume relation for h.c.p. Fe up to over 300 GPa 



36|,|9]. 

There is also very recent evidence for their accuracy in predicting the high-pressure phonon 

of DFT 



spectrum of h.c.p. Fe 143]. We use the ultra-soft pseudopotential implementation 
with plane-wave basis sets, an approach which has been demonstrated to give results for solid 



Fe that are virtually identical to those of all-electron DFT methods [39]. Our calculations are 
performed using the VASP code ]^5[, which is exceptionally stable and efficient for metals. We 
implemented a scheme for the extrapolation of the charge density which increases the efficiency 
of the molecular dynamics simulations by nearly a factor of two ]p6| . The technical details of 



pseudopotentials, plane- wave cut-offs, etc. are the same as in our previous work 



5.2 The hquid 

Our ab initio m.d. simulations on the liquid, which we used to calculate W{N, Nx) and hence 
the chemical potentials, were all performed on systems of 64 atoms, with a time-step of 1 fs 
and with F-point sampling of the electronic Brillouin zone. In our previous calculations on 
pure liquid Fe ]|l^, we showed that F-point sampling on a 67-atom cell underestimates the free 
energy by ca. 10 meV/atom; this is a completely negligible error for present purposes. The 
calculations were done at T = 7000 K and at the volume/atom V/N = 6.97 A/atom, which for 
pure Fe gives a pressure of 370 GPa. This pressure is somewhat higher than the ICB pressure 
of 330 GPa ]30]. The temperature is also higher than that at the ICB: our ab initio melting 
curve gives a melting temperature of ~ 6350 K (or ~ 6200 K after the correction due to our 



estimate of likely DFT errors 
than some other estimates 



15] at the ICB pressure of 330 GPa, which is already higher 
17]. But we shall see below that depression of freezing point 
due to impurity partitioning lowers this by a further 700 K. We have made rough estimates 
which show that the difference between 7000 K and our estimated ICB temperature is unlikely 
to change the chemical potentials of S and Si by more than 0.1 eV and that of O by more than 
0.3 eV, which will have no significant effect on our conclusions. The difference of pressures 
should also make little difference. 



We have used thermodynamic integration to calculate W{N,Nx) for the three solute el- 
ements S, Si and O for A'x = 3, 6 and 12, corresponding to mole fractions of 4.7, 9.4 and 
18.8 %. In doing this, we have aimed to choose the number of A values large enough and the 
duration of the simulation at each A value long enough to give a precision on W{N, Nx) /Nx of 
ca. 0.05 eV for S and Si and ca. 0.1 eV for O. To illustrate how the thermal average {Ui — Uq)x 
in Eq. ( p3D depends on A, we display this quantity in Fig. |l] at five equally spaced A values 
for the oxygen system with A^x = 12- We see that the dependence on A is not far from linear. 
Using Simpson's rule to perform the integral, we compared results for W{N, Nx)/Nx using the 
five A values 0.0, 0.25, 0.50, 0.75 and 1.00 with those obtained using only the three values 0.0, 
0.5 and 1.0, and found that they differ by less than the statistical error. Since the replacement 
of Fe by O is a greater perturbation than that of Fe by S or Si (see below), we have taken 
this as justification for using only three A values in all the thermodynamic integrations. Our 
numerical results for W{N,Nx)/Nx for X = S, Si and O, together with the linear least-square 
fit of Eq. (p^), are reported in Fig. ^, and the resulting values of a = fJ-xA ^^'^ ^ given in 
Table [l|. 



As explained in Sec. 4T, the b values have to be corrected as in Eq. (pq) in order to obtain 
Ax, and this requires the partial molar volumes vx- We obtain the partial molar volumes from 
the simulations just described by studying the pressure change resulting from the replacement 
of A^x atoms in the pure liquid by atoms of X at constant volume - this is straightforward, since 
the pressure is automatically calculated during the constant-volume simulations. We find that 
within the statistical errors the change of pressure is linear in cx for all three impurity species. 
We then use the fact that vx — va = (v / BT){dp/dcx)T', for Bt we use the bulk modulus of the 
pure liquid, which we know from our previous work |18|. The calculated vx values are 6.65, 
6.65 and 4.25 for S, Si and O respectively, compared with the volume per atom in the pure 
liquid of 6.97 A'^. We note that S and Si have almost exactly the same volume as Fe, but that 
the volume of O is considerably smaller. Using these vx values in Eq. (28), we now obtain the 
results for Ax given in Table [l|. We see that the difference between 26 and Ax is very small for 
sulfur and silicon, as expected, but is substantial for oxygen. 



5.3 The solid 

The available evidence strongly indicates that the stable crystal structure of Fe at the pressures 
and temperatures of the Earth's core is hexagonal close packed (h.c.p.) |4|], and this is the 
structure adopted in our calculations. We first present our calculations for S and Si, and then 
summarize briefly our results for the more complex case of O, which have already been reported 
elsewhere. 



5.3.1 Sulfur and silicon 

The calculations are performed on a 4x4x2 h.c.p. supercell containing 64 atoms, with a 3x3x2 
Monkhorst-Pack [^] grid of electronic /c-points which give free energies converged within a few 
meV/atom. In our calculations on the static zero-temperature lattice, we find that when a 
single Fe is replaced by S or Si at constant volume, the relaxation of neighboring atoms is very 
small, and the pressure change is also small. For the atomic volume v = 6.97 A^/atom, the 
partial molar volumes calculated without lattice vibrations give the differences vs — vpe = —0.32 
and vsi — fpe = —0.32 A'^, which are extremely small compared with upe- We assume that that 
these differences will not be significantly affected by thermal effects. 

We now turn to the harmonic frequencies uJn and a;^ needed for the harmonic difference of 
chemical potentials /^xa™ (^^^ (Po|))- We calculate these using a supercell of 64 atoms in all 
calculations. For the pure Fe system only two independent displacements of a single atom are 
needed to obtain the full force constant matrix We displace the atom by ca. 0.015 A in each 
direction, which is known to be small enough to ensure accurate linearity between forces and 
displacements. For the calculations where one Fe is substituted with S or Si, the symmetry 



is much reduced, and the number of atoms to be displaced is 15, with the total number of 
independent displacements being 33. For the systems with two solute atoms, the symmetry of 
the system is reduced even further, and we need to displace 20 atoms in all possible directions, 
for a total of 60 displacements. There are two distinct ways of putting two S or Si atoms on 
nearest-neighbor sites: the first has both sites in the same basal plane, and the second has 
them in adjacent basal places. Within our errors, we cannot detect the free energy difference 
between the two arrangements. Since the zero temperature value of the difference fx — fpe 
is very small, we have not attempted to calculate its high temperature value in the harmonic 
approximation, and we report in Table || the zero temperature value. The correction to Ax is 
negligible anyway and can be ignored. 

For sulfur and silicon we neglect anharmonic corrections. In our previous work on pure 
Fe Q we showed that at ICB conditions the anharmonic contribution to the free energy is 
roughly 60 meV/atom. In this case we are concerned with free energy differences between the 
pure Fe system and a system where one of the Fe atoms has been substituted with X, so the 
difference of the relative anharmonic contributions to the free energies is presumably smaller 
than that. 

Our calculated values of /x^a ^^'^ at w = 6.97 A^/atom and T = 7000 K are reported 
m Table 



5.3.2 Oxygen 



As emphasized above and in previous work |31|, substitutional O in h.c.p. Fe is highly anhar- 
monic, because O is considerably smaller than Fe and has great freedom of movement, so that 

ts 

the harmonic approximation is completely inadequate for calculating /^opg- We gave a brief 
summary in Sec. 4^ of the thermodynamic integration techniques used to do the calculations. 
The numerical result for fJ-Qp^ at v = 6.97 A^/atom and T = 7000 K is reported in Table |l|. 
We have not attempted to calculate Ax for X = O, since this would be extremely demanding, 
and turns out to be unnecessary in our analysis of core composition. 

To calculate vq — vpe in the solid we have repeated the calculations at different volumes 
and numerically differentiated the results, as described in section 4.3. The value of vq — vpe is 
reported in Table |l[ 



5.4 Core composition and temperature 

Some crucial features of our results are immediately clear from Table |l]: the liquid-solid dif- 
ference /ijf^ — /ixA is negative in all cases; its magnitude is somewhat smaller than k-sT for S 
and Si, but is much bigger than /cbT for O. This implies that the solutes will all partition from 
solid into liquid, as expected; but the partitioning will be weak for S and Si and very strong 
for O. 

To see the implications in more detail, consider the case of Fe/S. If we postulate that the 
core is an Fe/S binary alloy, then we can estimate the mole fraction Cg in the outer core by 
noting that the density of pure liquid iron at the ICB pressure is ca. 6 % higher than the 
values obtained from seismic data . We therefore add sulfur to the liquid until the density 
is reduced to the required value, which gives Cg = 0.16. Now if we ignored the dependence of 
/is on concentration, then the value //g — /ig* = —0.25 eV would give Cg/cg = 0.66, so that 
Cg = 0.11. However, the positive As values mean that both /ig and /ig increase strongly with 
increasing mole fraction of S, and this will tend to equalize the mole fractions in solid and liquid. 
If we solve Eq. (|^ self-consistently for Cg with the given Cg, we find Cg = 0.14. But a 14 % 
mole fraction of S in the inner core is completely incompatible with the seismic measurements. 
We can use the Cg volume together with the partial molar volume Vg to calculate the change 
of density of the solid due to dissolved S, and hence the ICB density discontinuity. We find 
the discontinuity is increased from the pure-Fe value of 1.8 % up to 2.7 it 0.5 %, which is still 



much less than the seismic value of 4.5 it 0.5 %. This means that the binary Fe/S alloy can 
be ruled out as a model for core composition. The argument is still stronger for Si, since the 
chemical potentials in solid and liquid are even more similar than for S. We conclude that the 
binary Fe/Si model must also be ruled out. 

For O, the situation is the opposite. The difference of chemical potentials in liquid and 
solid has the very large value /^Qp^ — ^Qp^ = —2.6 eV, which imples a strong partitioning 
from solid to liquid. If we repeat our analysis of the outer-core density with the partial molar 
volume Vq, we find that an oxygen mole fraction Cq = 0.18 is needed to match the density 
of the outer core. Eq. @ then gives Cq ~ 0.003, so that the O concentration in the inner 
core is very small. With our calculated Vq value, we then find an ICB density discontinuity of 
7.8 lb 0.2 %, which is markedly larger than the seismic value. A binary Fe/0 model can thus 
also be ruled out. 

Although all the binary models fail, the seismic data can clearly be accounted for by ternary 
or quaternary alloys of the three impurities. Ab initio calculations on such liquid and solid 
alloys would certainly be feasible with the methods we have developed, but would need a 
considerably greater effort. If we assume for the moment that the different impurities do not 
affect each other's chemical potentials, we can use our present results to construct a model for 
the core composition. We have seen that S and Si alone cannot explain the density jump at 
ICB, so there must be some O in the outer core. If we dissolve some O in liquid iron, together 
with S/Si, maintaining the density of the alloy equal to the density of the core, we increase 
the density jump at the ICB. This is because hardly any O goes into the solid. We therefore 
continue to add O till we match the density jump at ICB. The resulting chemical compositions 
of the inner and outer core are summarized in Table ^. 

With these compositions, Eq. (Q) now allows us to determine the shift of melting tem- 
perature from that of pure Fe at the ICB pressure; we find ATm = —700 it 100 K. Compar- 
ing this with our earlier ab initio melting temperature Tm = 6200 — 6350 K for pure Fe at 
p = 330 GPa |18], we obtain the estimate Tjcb ~ 5600 K for the temperature at the bound- 
ary between inner and outer core. This is quite close to estimates that have been made in 
other ways [^] . The implications of our temperature and chemical composition results for the 
understanding of the Earth's dynamics and past history will be explored elsewhere. 



6 Discussion and conclusions 

We have shown the practical feasibility of calculating completely ab initio chemical potentials 
in liquids and solids, and hence the ab initio treatment of chemical equilibrium between co- 
existing phases. The practical benefits of being able to do such calculations have also been 
illustrated by showing how they can help to improve our understanding of a controversial 
and important chemical-equilibrium problem in the earth sciences. We note that, although 
the calculations are demanding at present because of the need to perform substantial ab ini- 
tio molecular dynamics simulations, the underlying concepts are rather straightforward, and 
represent a simple extension of well-known classical techniques. 

In conclusion, we want to stress that the techniques should have rather wide applications. 
Although we have chosen to focus on the partitioning of impurities between coexisting solid 
and liquid phases, the methods could equally well be used to study partitioning between liquid 
phases, or between solid phases. The ability to calculate ab initio chemical potentials in liquids 
also makes it possible to contemplate the ab initio calculation of the solubility of solids, liquids 
or gases in liquids. The practical application of these ideas is likely to be limited only by the 
need to find economical thermodynamic integration paths for transforming chemical species 
into each other. 
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Appendix: from constant volume to constant pressure 

We explained in the text that the isobaric dependence of solute chemical potential on concen- 
tration is obtained from {dm/dcx)p, where m is the non-trivial part of the solute chemical 
potential, defined in Eq. (IT^). However, the quantity given by our ab initio calculations is 
{dm/dcx)v- We derive here the relation between the isobaric and isochoric derivatives of m. 
We start by noting that: 

{dm/dcx)p = {dm/dcx)v + {dm/dV)cy^{dV/dcx)p 

= {dm/dcx)v + {dm/dV)c^N{vx-VA) , (Al) 

with T held constant throughout, where is the total number of atoms, and we have used the 
basic definition of the partial molar volumes vx and f a of solute and solvent. Next, we refer 
to Eq. (15) to see that: 

{dm/dV)c^ = {difix - f^A)/dV)c^ , (A2) 

which can be reexpressed as: 

{dm/dV)c^ = (5(/xx - fXA)/dp)cAdp/dV)c^ = -{vx - va)Bt/V , (A3) 

with Bt the isothermal bulk modulus, and we have used the relations {dfj,x/dp)cx = vx and 
{d^A/dp)cj^ = VA- Combining Eqs. (Al) and ( |A3D , we have: 

{dm/dcx)p = {dm/dcx)v - nBrivx - va^ , (A4) 



where n = N/V is the overall atomic number density. 
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-0.32 


-0.32 


-2.35 




5.9 


2.7 




A^XA A^XA 


-0.25 ±0.04 


-0.05 ±0.02 


-2.6 ±0.2 



Table 1: Calculated chemical potentials (eV units) and partial atomic volumes fx (A^ units) 
of solutes X = S, Si and O in liquid and h.c.p. solid Fe at conditions close to those of the 
Earth's core (see text). Chemical potential of X is represented at low mole fraction cx by 
/^x = kBTlncx + fj-x, with fix linearized as /ix — /^x + ^xcx- The quantity ^j^^ is /jy^ — /Xpg, 
with //pg the chemical potential of pure solvent Fe; vxA is vx — vpe, with vpe the volume per 
atom in pure Fe. The meaning of the calculated quantity 6x used to obtain Ax is explained in 
Sec. [4.1| . Superscripts / and s indicate liquid and solid. 



Solid Liquid 
Sulfur /Silicon 8.5 ± 2.5 10 ± 2.5 
Oxygen 0.2 ± 0.1 8.0 ± 2.5 



Table 2: Estimated molar percentages of sulfur, silicon and oxygen in the Earth's solid inner 
core and liquid outer core obtained by combining ab initio calculations and seismic data. 
Sulfur /silicon entries refer to total percentages of sulfur and/or silicon. 




Figure 1: The integrand {Ui—Uq)\ (eV units) appearing in the thermodynamic integration used 
to calculate the free energy change W{N^ Nx) when A^x atoms of pure solvent are converted 
into solute atoms, with total number of atoms in the system = N (see Eq. (p3|)). Results 
shown refer to oxygen solute for A'x = 12 and N = 64. Filled circles show values computed 
from ab initio m.d. simulations, with bars indicating statistical errors. Curve is a polynomial 
fit to the computed values. 
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Figure 2: The calculated free energy change W{N, Nx) when Nx atoms of pure solvent arc 
converted into solute atoms, with total number of atoms in the system = N. Quantity plotted 
is W{N, Nx)/Nx (eV units) as a function of concentration cx = Nx/N for liquid and solid 
solutions of S, Si and O in Fe. Filled circles are results for liquid, with bars indicating statistical 
errors, and straight dotted line being a least-squares fit to these data. Continuous curves for 
S and Si show results for solid solution obtained from Monte Carlo calculations based on ab 
initio free energy of nearest-neighbor interaction. Open circle with error bar is result for O in 
solid from thermodynamic integration. 




